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Refraction of a Longuet-Higgins Gaussian sea by random ocean currents creates persistent local 
variations in average energy and wave action. These variations take the form of lumps or streaks, and 
they explicitly survive dispersion over wavelength and incoming wave propagation direction. Thus, 
the uniform sampling assumed in the venerable Longuet-Higgins theory does not apply following 
refraction by random currents. Proper handling of the non-uniform sampling results in greatly 
increased probability of freak wave formation. The present theory represents a synthesis of Longuet- 
Higgins Gaussian seas and the refraction model of White and Fornberg, which considered the effect 

OO of currents on a plane wave incident seaway. Using the linearized equations for deep ocean waves, 

we obtain quantitative predictions for the increased probability of freak wave formation when the 
refractive effects are taken into account. The crest height or wave height distribution depends 
primarily on the "freak index", 7, which measures the strength of refraction relative to the angular 
(-H spread of the incoming sea. Dramatic effects are obtained in the tail of this distribution even for 

the modest values of the freak index that are expected to occur commonly in nature. Extensive 

'""5 comparisons are made between the analytical description and numerical simulations. 

I 1 I. INTRODUCTION 

Q 

Freak waves in the ocean, known also as rogue waves or giant waves, are waves of extreme height relative to the 
^ typical wave in a given sea state, and may arise during a storm or in relatively calm seas. Due to their steepness, 
• ^ such waves pose great risk to cargo ships and even to large cruise liners. Well-publicized and documented encounters 
^ include a 25.6 meter wave that hit the Draupner oil platform in the North Sea in 1995, two ships that suffered damage 
' at 30 meters above sea level from a single wave in the South Atlantic in 2001, and the cruise liner Norwegian Dawn 
^ that met a series of three 21 meter waves off the coast of Georgia in 2005. In March 2007, the MS Prinsendam was 
hit by a 21 meter tall freak wave in the Antarctic. Satellite images taken over three weeks in 2001 and analyzed as 
part of the European Union Max Wave project [1] detected ten waves of height above 25 meters, suggesting that such 
T-H waves commonly occur in the world's oceans. 

Random (constructive) linear superposition of many plane waves with differing direction and wavelength, as in the 
Longuet-Higgins random seas model 0, offers a simple statistical explanation for the occurrence of freak waves. By 
^—i the central limit theorem, the sea surface height in this model must be a Gaussian random variable with some standard 
deviation cr. In the limit of a narrow frequency spectrum the crest height then follows a Rayleigh distribution: the 
g probability of crest height exceeding is given by 

> -PRaylcigh(i/) = e-^'/2-' . (1) 
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Note that for linear waves there is an exact symmetry between crests and troughs, so that for a crest height of H, the 
^ wave height (crest to trough) is given by 2H. Conventionally, a freak wave is defined as H > AAa, or 2H > 2.2 SWH, 
where the significant wave height SWH w 4. Oct is the average of the largest one third of wave heights in a time 
series [3]. The Rayleigh distribution thus predicts freak waves to occur with probability 6.3 • 10~^ and extreme freak 
waves of crest height H > 6a (or 2H > 3 SWH) to occur with probability 1.5 • 10^*. Observational data [T] suggests 
that this purely stochastic Rayleigh model significantly underestimates the actual number of freak waves. A review 
of several alternative theories of the freak wave phenomenon appears in Ref . [4j . 

Nonlinear instability effects |51 have been extensively and successfully studied as a mechanism of freak wave 
formation. However, the effects of such instabilities depend sensitively on initial conditions, and the full numerical 
computations starting from a generic random sea state are costly. Thus, it is difficult to obtain quantitative predictions 
of the crest height distribution, analogous to ([T]), except in approximations such as the Nonlinear Schrodinger Equation 
(NLS) or the Dysthe equation [SI |7], which are valid for small to moderate values of the wave steepness. Since 
nonlinear effects scale as powers of the wave steepness kH, where k is the wave number, strongly nonlinear evolution 
is more likely to be triggered in an initial condition where the waves are already unusually high. Thus the tail of the 
crest height distribution may be dominated by linear triggering mechanisms that produce slightly taller-than-average 
waves especially conducive to nonlinear instability. Identifying such triggering mechanisms, and combining them with 
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nonlinear evolution, would allow for quantitative predictions of the freak wave probability distribution, and would 
also advance the long-term goal of freak wave forecasting. 

One such linear triggering mechanism for deep water waves is the focusing or refraction of an incoming plane wave 
by random current eddies, as discussed by several authors including Peregrine [8, and White and Fornberg [9J, and 
motivated by the fact that many freak waves have been observed in regions of strong ocean currents. This mechanism 
must contain part of the story, since ocean currents do refract waves in the ocean. The difficulty is with the assumption 
of a single plane wave incident on the refracting currents. A plane wave leads to caustics or singularities with infinite 
ray density (smoothed out only at the wavelength scale), and consequently a repeated and reproducible pattern of 
freak waves, which will always appear whenever focusing currents are present. Statistical predictions have no place in 
such a theory. Of course, an incoming plane wave is a physically unrealistic model of most sea states, since directional 
and wave number spread in the incoming sea will smear out the above-mentioned singularities. 

These difficulties may be overcome by combining the stochastic random seas approach and the focusing approach, 
to obtain the distribution of crest heights resulting from a random incoming sea incident on a region of random eddy 
currents. Averaging over a random incoming sea smears out but does not entirely destroy the effects of the refraction. 
Lumps in wave action persist, varying from place to place by factors of 2 to 4 typically. For realistic values of the 
parameters, the predicted probability distribution of crest heights depends on a single quantity, the "freak index," 
which is a simple function of mean wave speed, rms current speed, and angular spread of the incoming waves. We shall 
see that realistic sea states can lead to enhancements of a factor of 50 or more in the probability of forming freak waves, 
over and above purely Gaussian seas with the same averaged action density. The greatest probability enhancement 
occurs for the largest freak waves, while the probabilities associated with typical wave heights are scarcely affected. 

In Sec. [njwe review the focusing model, combine it with a stochastic incoming sea, and obtain analytical results 
for the probability distribution of crest heights when the freak index is small. The predictions are tested in Sec. |III| 
using numerical simulations in the ray limit (where the wavelength is small compared to the size of the eddies). The 
quantitative correspondence between ray dynamics and wave dynamics is discussed further in Sec. |IV[ in connection 
with analogous correspondence for the Schrodinger equation. We emphasize that the probability distributions obtained 
in Sec. [IT] and verified numerically in Sec. |III| are not the final word, but must rather be used as input to the full 
nonlinear evolution. This and other outstanding questions are explored in Sec. |VI[ 



II. THEORY 
A. Refraction 



We begin by briefly reviewing the focusing of a plane wave by a random current field. A detailed discussion may 
be found in the work of White and Fornberg [3] . 

The ray dynamics of deep-water surface gravity waves is governed by the usual eikonal equations for ray position 
r and wave vector k, 

dk duj dr duo 

dt ^ dt^ dji' ^ ' 

Here the dispersion relation is given by 

uj{k,r)^^/^\ + k-U{r), (3) 

where U(f) is the time-independent current velocity, assumed to be slowly varying on the scale of a wavelength. In 
the following analysis, we consider J7(r) to be a random field, with zero average velocity, typical velocity fluctuations 
of size uq, and spatial correlations on a distance scale ^. In the ocean, the typical eddy size ^ may vary between 20 
and 100 km, while the typical speed uq is generally below 1 m/s. 

An incoming plane wave of frequency lu moving in the y direction may be represented in phase space by initial 
conditions r = {x, 0) for all x and k = (0, k), where k = 27r/A is related as usual to the wave speed hy v = doj/dk — 
(.g/4fc)i/2 and to the wave period as T = 27r/w = 27r(gfc)i/2 (here we neglect the current velocity U , which generally 
is much less than w, but the exact expression ^ is used to relate lo, v, and k in numerical simulations, causing v and 
k to be spatially varying even for a constant frequency co). A typical wave period for deep-water ocean waves is 10 s, 
corresponding to wavelength A = 156 m and wave speed v = 7.81 m/s. 

When this incoming plane wave impinges on a random current field, it will undergo small-angle scattering, with 
scattering angle ~ mq/v after traveling one correlation length ^ in the forward direction. Eventually, singularities 
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appear that are characterized in the surface of section map [a:(0), fc2;(0)] — > [x{y),kx{y)] by 5x{y) /Sx{0) = 0, i.e., by 
local focusing of the manifold of initial conditions at a point. The formation mechanism of these singularities may 
be ascribed to a 'bad lens.' Whereas a good lens without aberration focuses all parallel incoming rays to one point, 
a bad lens only focuses at each point an infinitesimal neighborhood of nearby rays, so that different neighborhoods 
get focused at different places as the phase-space manifold evolves forward in y, resulting in lines, or branches, of 
singularities. The typical pattern is an isolated cusp singularity, S'^x{y) /Sx{0)'^ — 0, followed by two branches of fold 
singularities, as shown in Fig. [T] 




FIG. 1: A cusp singularity, followed by two branches of fold singularities, is formed as initially parallel rays pass through a 
focusing region. The two branches appear because the focal distance varies with the distance of approach from the center, as 
in a 'bad' lens with strong spherical aberration. After averaging over incident directions, the singularities will be softened but 
not washed away completely. 

The first cusp singularities appear after a median travel distance y — L ^ (.{uo/v)^^^^ 3> through the random 
eddy field [9l [10], when the typical ray excursion in the transverse x direction becomes of order ^. For realistic 
parameters, L 100 km or more is typical. [Note that the weakness of the scattering requires a wave to pass 
over many uncorrelated eddies before the first singularities are formed. This allows ray propagation to be described 
statistically without regard to detailed structure of individual eddies, using only the length scale ^ and dimensionless 
velocity ratio uq/v.] The singularities formed in this way are separated by distance ^ ^ in the transverse direction, 
and the typical deflection angle by the time these singularities appear scales as 

50«^^K/«)2/3. (4) 

The quantity 69 may be defined precisely as the rms value of the deflection angle evaluated at the forward distance 
L, where the averaging is performed over all rays and over an ensemble of random eddy flelds. The typical deflection 
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angle 56 does not depend on the size of the eddies but only on the velocity ratio uq/v: faster currents cause larger 
deflection. For example, for v = 7.81 m/s and uq = 0.5 m/s, we find (uq/v)'^^^ = 0.16, SO = 18°, and the median 
distance to the first singularity is i = 7.5^ (150 km if we take the eddy correlation length ^ to be 20 km). 

Of course, the singularities are a mathematical construct existing only in the ray limit k^^ — > oo. For wave 
dynamics, any such singularities must be smeared out at least on the scale of a wavelength, or more precisely, on the 
scale fc^T^ ~ k~^{uo/v)~^^'^ in the transverse direction. For example, a fold singularity will be softened over a distance 
scale k~^{kxS,y^^ [TT], resulting in a maximum wave intensity enhanced by a finite factor (k^S^)^^^ ^ (fc^)^/'^(Mo/w)^/^ 
compared to the background intensity, which is a factor of 5 or more for reasonable-sized eddies. Although infinities 
are thus absent from the wave dynamics, this model is still unsatisfactory, as it predicts that extreme freak waves will 
always occur whenever an incoming plane wave encounters a random eddy field. 

Subsequent propagation through the random current field produces new singularities, a nd th e resulting number of 



branches grows as e^^^ while the wave continues to travel forward. As we will see in Sec. II B below, wave intensity 
peaks associated with these later singularities become increasingly washed out due to growing spread in wave direction. 
Thus, the early generations of caustics, appearing soon after the incoming wave is first scattered by the eddy currents, 
will typically dominate the freak wave distribution. 



B. Averaging 

We now consider the more realistic situation where the initial plane wave is replaced by a random superposition of 
waves with a finite angular spread A9. Typically, Ad may take values of 10 to 30 degrees. Again, for k^S, ^ 1, we 
are justified in modeling the wave dynamics using the ray approximation, where the initial positions are still given 
^ — (2^,0) uniformly distributed over all x and the initial wave vectors are given hy k — (k sin , k cos 9) , with 
a Gaussian directional spread P{0) ~ /2(Ae) ^ pj-^ gg^,^ jy confirm the quantitative correspondence between 



wave and ray intensity statistics in the context of the linear Schrodinger equation.] 

This set of initial rays is refracted in its evolution through the eddy field, and by the time the first singularities 
appear, the typical ray is scattered by a transverse wave vector of size Sk^, as described by Q. Due to the finite 
spread Ak^ ~ kA6 of initial conditions, the singularities in position space are smoothed out. Instead we obtain 
regions of above average intensity near the positions of the would-be singularities, and regions of below average 
intensity elsewhere, as indicated in Fig. |2] The typical contrast between high and low intensity regions is determined 
by the dimensionless freak index [T^ . 

^ A0^ Ak^^ A0 ' ^ ' 

For example, when 7 <C 1 , scattering produces small fiuctuations of order Sk^ in the vertical thickness of the phase 
space region representing the evolved sea state in Fig. |2](Left). The mean vertical thickness of this phase space region 
remains Ak^, and the contrast ratio between maximum and average intensity therefore scales as 1 -I- 0(7). In the 
opposite limit, 7^1, the singularities in the ray density are only slightly softened by the spread Ak^ in the initial 
sea state, and the maximum intensity is enhanced by a factor scaling as 7"^/^ compared to the average intensity [13,. 

Thus, the likelihood of freak wave formation is greatest for small initial directional spread Akx and large defiections 
6kx, i.e., for a well-coUimated (long-crested) sea encountering a very strong random current field. On the other hand, 
we will see below that typical parameter values occurring in nature correspond rather to small or moderate values 
of the freak index, 7 < 1 or 7 ~ 1, where the small-7 expansion will be more relevant. Interestingly, even in this 
regime where most refraction effects have been washed out, one can nevertheless observe dramatic consequences for 
the probability of freak wave formation, particularly for the extreme freak waves. 

What justifies our focus on intensity fiuctuations associated with the first generation of singularities, those formed 
after typical travel distance L? Subsequent evolution through the random eddy field causes the direction of travel to 
undergo diffusion, so that the initial angular spread A9 becomes 

(A%))2«(A^)2 + O((50)Vi) (6) 
after distance y, and thus the effective freak index decreases as 

(7) 

at large distances y. The first generation of singularities, formed soon after entering a region of strong random 
currents, is responsible for producing the largest intensity fiuctuations and thus appears to be the most hospitable 
environment for freak wave formation. 
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FIG. 2: Left: An initial sea state characterized by a spread Ak^ ~ k/S.9 in the transverse wave vector is described by a strip 
of finite thickness in the {x,kx) phase space a.t y = Q (top). Under subsequent refraction through the eddy field, the sea 
state acquires additional fluctuations of typical size 5kx ~ kS6 as it evolves forward in y, and eventually singularities in the 
ray dynamics are encountered. Although all singularities are washed out due to the finite initial spread /S.kx, a pattern of 
high intensity (focusing) regions and low intensity (defocusing) regions is evident when the evolved phase space distribution is 
projected onto position space. Right: Further refraction through the random eddy field produces multiple tendrils in the phase 
space distribution and many high and low intensity regions in position space, which may be described statistically. 

We now consider the implications of these spatial variations in ray density for the distribution of crest heights. 
Since the ray density fluctuations occur on a scale larger than a wavelength, the local wave dynamics may still be 
described by a Gaussian-random Longuet-Higgins model, and M generalizes to 



Here d''^{x,y) — cr^/iay (a^, y), where Ij:g,y{x,y) is the local density of rays, normalized to unity for the initial sea state 
before scattering, and is the elevation variance associated with the initial sea state. 

Here we should note that both wave energy density and wave action density are proportional to the mean squared 
wave amplitude, which in turn corresponds to ray density. Thus the lumps we find in ray path density or squared 
wave amplitude imply, and can be converted into, lumps in either energy density or in action density. The conserved 
quantity in the presence of currents is wave action, rather than wave energy jl4j . In the present circumstances, the 
conversion factor relating action and energy is nearly constant: the effect of currents of strength 0.5 m/sec as against 
a typical 10 m/sec wave velocity gives corrections on the order of 5%. 

Now let TZa be the ratio of the local probability oi an H = aa event to the probability of the same event using a 
Rayleigh distribution. Then 



Eq. ([9| shows that local enhancement and suppression of the freak wave formation probability is very significant: in 
a zone where the local energy is just 50% above the mean (/ = 1.5), the frequency of a 4.4ct wave crest (the threshold 
for a rogue wave) increases by a factor of 25, while the frequency of a 6tT wave crest is enhanced by a factor of 400. 
At the same time, the low energy zones are remarkably quiescent: 4.4(t events are 20 times less likely in a patch down 
only 25% in energy density from the mean (/ = 0.75). 



P\oca\{H) 



e 



(8) 




(9) 
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After spatial averaging, we obtain the total probability of a crest height exceeding iJ, 

P(H) = J dIV{I) e-^V2^^/ ^ (10) 

where V{I) ^ J dxdy S{I — /ray (2;, y)) is the probability distribution of energy densities obtained using ray dynamics. 
This local rescaling of the random wave model is unproblematic for linear wave equations, and has been used suc- 



cessfully in computing wave function statistics for the linear Schrodinger equation ^\ [TS] (see also Sec. IV ) . Caution 

must be used in applying such rescaling to a nonlinear wave equation, as we know that over sufficiently long time 
scales, nonlinear instabilities will cause the sea state to re-equilibrate to a longer or shorter mean wavelength after a 
change in the energy density. This and related issues are addressed in Sec. |VI[ 



C. Analytical Results 

We work in the regime of small or moderate freak index 7, where scattering effects are relatively weak compared 
with the angular spread of the incoming waves. This limit is most likely to be encountered in nature, and yet yields 
surprisingly large effects in the tail of the freak wave probability distribution. Consider a Gaussian-distributed local 
energy density, with a standard deviation e proportional to 7: 

-(/-1)V2£= 



where 



Then the distribution of crest heights becomes 



6 = 07. (12) 



g-(/-l)V2£%-ffV2o-^/ 



P(H)^J4I . (13) 

For small fluctuations e, or large heights H, the integral may be evaluated by the method of steepest descent, expanding 
around the maximum I — I ~ z, and we obtain 

PiH) = ^li^/^'^'^''^'^^'^ > (14) 
where z is defined implicitly as a function of H through the cubic equation 



z{l + zf^^-^. (15) 
The validity of the steepest descent approximation requires either small density fluctuations e <C 1 or very tall 



waves H/a ^ 1, or both, but the right hand side of (15 1 may in general be of order unity. However, if we consider 
the probability of waves of a fixed height H in the limiting case of a very small freak index, e — > 0, then z becomes 
small and we obtain the perturbative result, 



Ppcrt(i/) = 



l + e--V-%0(e^) (16) 



i.e., a polynomial multiplying the original Rayleigh distribution. This perturbative result is analogous to quantum 
wave function intensity distributions in the presence of weak disorder or weak scarring by periodic orbits |16L I17j . 
According to the perturbative expression, the probability of seeing a wave height equal to or greater than 2H = Aa 
(the significant wave height of the initial sea state), is unaffected by refraction, but the probability is enhanced as we 
consider more and more extreme waves. We note that PpertiH) is a cumulative distribution function. The probability 
density of encountering a wave crest of height H is given by dPpcTt{H) / dH , which is enhanced both for very small 
heights {H < 1.17 cr) and very large heights {H > 2.22 a) and is reduced for waves of intermediate height. This is 
consistent with the physics of energy focusing, which increases fluctuations while keeping the average energy density 
unchanged. 
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On the other hand, if we consider the extreme tail H/a — > oo for a given freak index, the right hand side of (15 1 
becomes large, and so does the z parameter. We obtain 



Pt&iAH) ^ exp 



2 V 2ct2 / V 2ct2 / '3 



(17) 



Note that the leading term in the exponent grows only as H^^^ for very large crest heights H, in contrast with the 
much faster growth for the Rayleigh distribution. However, depending on the value of the freak index 7 (or 
equivalently on the typical intensity fluctuation parameter e), the probabilities associated with this asymptotic tail 
may be too small to be observable in practice. 

Also of interest are the moments of the crest height distribution, ff^n Since the crest heights are locally Rayleigh- 
distributed, we may write 

H"^ ^ -ffRaylcigh-^ray(a;, v) (18) 

where ^^Rayicigh ^ Rayleigh-distributed random variable corresponding to the average surface elevation variance cr^ , 
and /ray (a;,?/) is a slowly changing scaling factor that describes energy density variations. Averaging over space. 



tr2n _ Ej2n . Jn 

~ -"Rayleigh ^ 



^ravMgh • dIV{I)I\ (19) 



Rayleigh 

where the first factor is the corresponding moment in the absence of refraction, while the second factor depends only 
on the ray dynamics and takes into account fluctuations in the local energy density. Assuming once again a Gaussian 



distribution of ray densities for small freak index (111, we have 



(/-l)» = e"(n-l)!! (20) 

for even n and for odd n, and in particular the scaling of the even moments with the freak index 7 is predicted to 
be 



i?«-((/-l)")*-/|-«/f- 



7. (21) 



In a regime where the Gaussian approximation is inadequate, the moments /" or i?" = (/ — 1)" may be evaluated 
directly from a numerical ray dynamics simulation. We note that / = 1 is required by probability conservation in 
the ray dynamics, while the higher moments /" for n>2 are necessarily greater than 1 in the presence of refraction. 
However, for small freak index 7 <C 1, the corrections to Rayleigh are tiny except for the very high moments, n ~ I/7 
for /" or n 1/7^ for i?„. Only in these high-order moments can one clearly see the dramatic refraction-induced 
effects which dominate the tail of the crest height distribution. 

Finally, we must consider the effect of a flnite range of wave frequencies (or speeds or wavelengths) in the initial 
sea state, as given for example by the JONSWAP spectrum [T^. For simplicity of presentation, consider flrst the 
simplifled case where the energy is uniformly distributed among all wave speeds in the interval [v — At;, v + Au]. The 
total finite-bandwidth ray intensity is then given by 



1 



Ifh{x,y) - — / dv'h,{x,y), (22) 
2Au /„_Ai, 

where ly'ix, y) is the ray density associated with the single wave speed v' , and obeys the statistical properties discussed 
above for freak index 7 ~ {uq/v'Y^^. We are of course interested in fluctuations of the total ray intensity around its 
average, e.g., the variance 



dv' I dv" Slv>{x,y)SIv"{x,y) , (23) 



^ (2Ai;)2 J^^Av 
where 

(x, 2/) = (x, y) - 1 ^ 7 ^ . (24) 
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Two effects must be taken into account in determining the effect of finite but small Av in (23 1. First, the size of 
the intensity fluctuations at a given velocity is velocity-dependent, i.e., SIy'{x,y)SIvi{x,y) ^ w'""*/^, and since the 
function v'~^^^ has a positive second derivative, integrating over a s ymm etric interval of size Av around the central 



velocity v will lead to an 0((Aw/w)^) enhancement in the variance (23) compared with the result for a single wave 
velocity v. Secondly, replacing the initial velocity v with v' causes the rays to deviate by a typical distance ~ ~") ^ 
by the time the first caustic is reached, and thus the intensity fluctuation map SIy>{x^y) is shifted with respect to 
the original map SIy{x,y) by a similar displacement. Now assuming that the original ray intensity map SIv{x,y) is 
sufficiently smooth on the eddy scale ^ (true for large angular spread A9 or small freak index 7, as discussed above), 
we find that the correlation between 5Iv{x,y) and 6Iy'{x,y) should fall off as 1 — 0{{v' — vY /v^). Both effects are 
second order in Av, and combining them wc find that the variance of the energy density for finite bandwidth behaves 
as 

el^e^[l±0{Av/vf], (25) 

where e is the zero-bandwidth result evaluated using the central velocity v, where the sign and numerical coefficient 
of the leading correction may be 7-dependent. Replacing the uniform velocity spread with a more realistic spectrum, 
such as JONS WAP, will merely change the numerical coefficient. In Sec. |III| wc will confirm that the quantitative 
size of the finite-bandwidth effect is very small in practice, and will almost certainly be overwhelmed by nonlinear 
corrections that are explicitly excluded from our model. Thus, within the linear approximation one is generally well 
justified in ignoring finite-bandwidth effects, and simply relying on the freak index obtained from the mean wave 
speed. 

III. NUMERICAL SIMULATIONS 

Following the work of White and Fornberg [9^ , for our numerical simulations we generate incompressible random 
current fields U{r) as 

U.,{f) = -d4,{T)/dy- Uy{f) = di^{r)/dx . (26) 
Here the two-dimensional stream function ijj{'r) is Gaussian distributed with Gaussian decay of spatial correlations: 

iAf) = 0; iA/)iJ{r^ ~ e-^^"'-^"')'/^?^ ^ (27) 



and normalized so that |[7(r)P — Uq. The choice of a Gaussian-distributed and Gaussian-correlated stream function 
is made for convenience and for consistency with Ref . \^ . The theoretical discussion in the previous sections does not 
depend on any specific choice of a random ensemble, but only on the length scale ^ and velocity scale uq. 

The calculation is performed on a 640 km by 640 km grid, with an eddy correlation length ^ = 20 km and periodic 
boundary conditions in the transverse {x) direction. Without loss of generality, the rms current speed uq is set to 
0.5 m/s. Again, the specific values of ^ and uq are arbitrary and serve merely to set the scale for the simulation. The 
refraction strength as measured by SO or Sk^ may be controlled by varying the incoming wave velocity v, in accordance 
with ([4]), while the angular spread AO is controlled directly in the initial conditions. To avoid boundary effects, ray 
trajectories are launched at a distance yo = 50 km inside the random eddy field, uniformly spaced in the transverse x 
direction, and with a range of wave vector directions Oq. The initial wave vector for each ray is fco = (A;o sin Oq, kg cos Oq), 
where the initial wave number fco is chosen to correspond to constant frequency uj = 27r/(10 sec) in accordance with 
These trajectories are interpolated and weighted with ^(6*0) ~ e~^o/^^^^^^ to produce a ray density map I^ayix, y). 
Typical density maps for the same random eddy field but two different values of the initial angular spread are shown 
in Fig. [3j 

Recalling that the median distance at which initially parallel rays form a caustic is approximately 150 km for our 
parameters, we collect statistics on ray density I,-s^y(x,y) over the region j/o + 12.5 km < y < yo + 250 km, which 
includes the most significant density ffuctuations. Several moments of the ray density distribution for various initial 
angular spreads AO are shown in Fig. |4] The initial wave velocity v and rms current speed uq are kept fixed in these 
simulations, so varying AO corresponds to a variation of the freak index 7 ([s]). 

We notice first that the norm R2 is simply the standard deviation e of the ray intensities, and its linear scaling 



with the freak index 7 for small 7 confirms the prediction ( 12 1. Similarly the scaling of the higher moments with 7 is 



consistent with the prediction of (20 1 and (21). The choice a = 0.25 for the proportionality constant in (12) results 



in good agreement with (21) for all moments in the regime of small or moderate freak index, as indicated by the 
theoretical lines for R2, R4, Rio, and i?20 in Fig- 14] The fact that all moments are described by a single constant 



9 




FIG. 3: A ray density map is shown for rays moving through a 640 km by 640 km random eddy field with rms current 
Mo = 0.5 m/s. The rays are coming in from the top with wave velocity v = 7.81 m/s, uniform initial distribution in the x 
direction, and initial angular spread A9. Top: = 5°, corresponding to a very high freak index 7 = 3.6. Bottom: = 25°, 
or 7 = 0.7. 
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FIG. 4: The norm 
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of the ray density I^a,y{x,y) (Fig. |3j) is calculated for several n and for different initial 
angular spreads A9 (corresponding to different values of the freak index 7). From left to right, the data points represent 
A6 = 25°, 20°, 15°, 10°, and 5°. In each case, the intensities are sampled on a uniform rectangular grid extending in the 
longitudinal y direction from a distance of 12.5 km to a distance of 250 km beyond the y— position where the rays are initially 
launched. The theoretical lines for R2, R4, Rio, and R20 are obtained from (121, (201. 
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FIG. 5: The fraction P{H) of crest heights greater H is computed from the ray densities /ray(a;, y) using the assumption 
of locally Rayleigh fluctuations (solid lines). From bottom to top, the five solid lines show results for initial angular spread 
AS = 25°, 20°, 15°, 10°, and 5°. The small-7 analytic prediction (141 is shown by a dashed line in each of the five cases. The 
Rayleigh distribution ([l]), which describes the refraction- free 7^0 limit, is also shown for comparison. 



a validates the assumption of Gaussian density fluctuations in this regime. Notably, linear dependence on the freak 
index 7, especially for the high moments that govern the tail of the crest height distribution, is evident for all but the 
largest 7 (corresponding to initial angular spread — 5°). Thus, the small-7 approximation appears to be justified 
for most sea conditions expected in nature. Furthermore, the uppermost line in Fig. |4] shows that the maximum 
intensity fluctuation R^o — max(|/ — 1|) also scales linearly with 7, as expected for a finite sample. 

Given the ray densities I^a,y{x, y), the probability P{H) of encountering acrest height larger than H may be obtained 
using (10 1. The results for several values of the initial angular spread A9 are shown in Fig. [s] and compared with 
the analytical prediction ( 14 1 , which was derived in the limit of very small freak index. The Rayleigh prediction 
([ij serves as a baseline for comparison. In Fig. |6] we calculate the factor by which P{H) is enhanced over the 
Rayleigh prediction, for several values of H . We see that under realistic conditions, refractive effects may enhance the 
probability of encountering a freak wave {H = AAa) by as much as an order of magnitude, while the probability of 
extreme freak waves {H = 6(t) may be enhanced by more than two orders of magnitude, depending on A9. The small- 
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FIG. 6: The increase in the number of waves with crest height greater than H is shown as a function of freak index 7. Solid 
squares correspond to H = AAa (the traditional definition of freak waves), solid circles correspond to H = 5a, and solid 
triangles correspond to H — 6a (extreme freak waves). The open symbols indicate an analogous calculation, but for finite 
bandwidth Av/v = 0.4. From left to right, the data points represent initial angular spread A6 = 25°, 20°, 15°, 10°, and 5°. 
The dashed lines are obtained from the small-7 analytic prediction (141. 



7 analytic approximation provides a reasonable estimate of this enhancement for all but the least realistic situation 
of very long-crested incoming waves {A9 = 5°, or 7 = 3.6). 
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FIG. 7: The standard deviation R2 = and maximum fiuctuation Roc = max |/ — 1| of the ray intensity / are plotted as 

m Fig. |4] but varying the velocity v of the incoming wave as well as the angular spread AQ of the incoming sea. The velocities 
chosen correspond to « = /3wo, where /3 = 0.8, 1.0, 1.2, and 1.4, and wo = 7.81 m/s is the original velocity used in the previous 
figures. The five symbols of each type correspond to the same five values of AQ as in Figs. [4] and |5] 

Until now, we have been varying the freak index 7 = 59 / AO by adjusting only the angular spread AO of the 
incoming sea. To verify that the freak index is in fact the single parameter determining the probability of freak wave 
formation in our model, we must likewise vary the refraction strength 50 ~ (uq/w)^/'^. Since the eikonal equations ([2]), 
([3]) are manifestly invariant under the simultaneous rescaling of the wave speed v and rms current speed Uq {v —>■ qv, 
uo quQ, k q^^k, oj q^^to), we may without loss of generality vary v while keeping the currents fixed. The 
results of such a calculation are shown in Fig. [7] and confirm that the energy density fluctuations leading to freak 
waves exhibit single-parameter scaling with the freak index 7. 

Of course, a realistic initial sea state has a finite frequency bandwidth, and consequently includes a range of 
wave speeds rather than a single speed v. Repeating our simulations for velocity uniformly distributed in the range 
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[v — Av,v + Av], we find little change in the size of ray density fluctuations, as compared with the monochromatic case 
of a single speed v. This is consistent with the discussion leading to (25 1, which showed that finite-bandwidth effects 
vanish at leading order in Av/v. For example, in the typical case of freak index 7 — 1.2 (central velocity v — 7.81 m/s, 
rms current speed uq = 0.5 m/s, and directional spread A9 = 15°), the moderate bandwidth Av/v — 0.2 results in 
an increase of less than 1% in the width efb of the intensity distribution, as compared with the monochromatic result; 
doubling the bandwidth to Av/v = 0.4 still produces only a 3% increase in the width of the intensity distribution. 
The effect on freak wave formation, for Av/v = 0.4, is indicated by the open symbols in Fig. [6] (the effect of bandwidth 
Av/v = 0.2 would be too small to be clearly visible on the scale of the figure). 

Another phenomenon observed in the numerical simulations, which may be of interest when comparing with obser- 
vational data, is the rapid change in mean wave direction that typically occurs in conjunction with hot spot formation. 
This is easy to understand in the large 7 limit, where a fold singularity results in a discontinuity in the ray density 
necessarily associated with a discontinuity in the mean ray direction (since the wave vector k associated with the 
fold will be different from the wave vector associated with other parts of the phase space manifold that project onto 
the same position f). For finite initial angular spread AO, such discontinuities are smoothed out as we have seen, 
but the residual effect remains. For example, on a square grid of cell size 1.25 km (« 8 wavelengths), we observe a 
maximum cell-to-cell change of 25° to 28° in the mean wave direction when 10° < A9 < 25°, to be compared with a 
median cell-to-cell fluctuation of 2° to 3° over a 640 km by 640 km grid for the same parameters. This rapid change in 
wave direction associated with entering a hot spot may be consistent with mariners' observation of some freak waves 
appearing at large angles from the mean wave direction. 



IV. SCHRODINGER EQUATION SIMULATIONS AND STATISTICS 

We now describe linear Schrodinger equation simulations, which are designed to check the assumptions and ideas 
presented so far: 1) Is the connection between mean squared wave amplitude and ray density sound, under the 
conditions of averaging over wavelength and especially direction? 2) is it correct to derive the global wave statistics 
as an integral over locally Gaussian seas (the local Rayleigh approach, see (|8|)? The linear Schrodinger equation with 
a random potential is a slightly imperfect laboratory for making comparisons with water waves, and of course the 
correct ray dynamics using eddy fields has already been presented above. However, solving the correct nonlinear water 
wave equations with proper dispersion on these eddy fields is difficult, and although we hope to do so in the future, we 
argue that the linear Schrodinger equation is adequate for the present statistical checks on ray-wave correspondence. 

In the wave simulations, statistics for rare events have to be gathered over large runs in space and time. One cannot 
sit only at the focal region of an eddy and have event after event; one has to wait and measure the statistics there as 
elsewhere. If we had used a plane wave, the wave would repeat itself periodically. With a random superposition of 
hundreds of waves, we never encounter repetitions in the time allotted to the simulations. 

The simulations are performed as follows. Both the ray and wave simulations use of course the same random 
potential, which plays the role of the random eddy field. [Note that the formation of caustics and energy lumps is 
independent of the microscopic mechanism of ray defiection; only the correlation length and mean deflection angle 
are important.] For each run, we produce a Gaussian random potential fleld V{r) by superposing 40 sine waves with 
random direction, phase, and amplitude. An overall scale factor controls the strength of the potential, and thus 
the distance to the flrst focal caustics. The random potential may be suppressed by a smooth tanh function in the 
entrance region (top) and/or exit region (bottom), to simulate eddy- free zones through which the random waves travel 
before and after refraction by the eddies. The images shown below are about 8 correlation lengths across. 

For the ray studies, we establish a rectangular position grid. Typically 40,000 ray trajectories are launched uniformly 
along a line in the zero-potential entrance region, with initial direction and speed taken from appropriate Gaussian 
distributions. To avoid pixelization effects, each trajectory is treated as a very narrow Gaussian density (typically 15 
pixels across, which is much less than the size of the energy lumps). 

On the same potential fleld we propagate waves through the region of study, using a 512 by 2048 square grid 
corresponding approximately to 80 by 320 wavelengths. Each incident waveset begins as an appropriate random 
superposition of 400-700 traveling plane waves, and is then propagated by the split operator fast Fourier transform 
method [19 (SOFT). Statistics are collected at every time step. The mean squared amplitude density plots and wave 
statistics data are obtained in each run (i.e., for a given potential, dispersion of incident directions, etc.) by repeating 
this process for 5 to 10 incident wavesets. 

Figures |8] and [9] reveal a great deal about the calculations and the results. Since our aim is to check the accuracy 
of freak wave statistics predicted from ray data, all runs are performed at moderate to large freak index 7, where 
statistically signiflcant numbers of freak wave events appear, and where ray-wave correspondence is most likely to be 
suspect. Each run takes a few hours on a single processor Macintosh G4 workstation. In Fig. [8] we see ray density 
(panels A, B) and wave intensity (C, D, E) data, together with the potential superimposed in E, for a freak index 
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7 « 3.4. Very impressive and detailed agreement is seen comparing the ray density and wave intensity averages; 
compare B (ray) with C (wave) especially. In A, the singular ray density for = (7 = 00, corresponding to a 
single incident plane wave, i.e., the White and Fornberg limit [S]) is shown in red for comparison, superimposed on 
the direction-averaged ray density for 7 « 3.4. In panel D, freak events of 4.4ct are superimposed in blue, and 6ct 
events (truly disastrous!) are shown in red. 
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FIG. 8: Ray density (panels A, B) and wave intensity (C, D, E) data, together with the potential superimposed in E, and 
freak events superimposed in D, for a freak index 7 « 3.4. Rays and waves are launched from the top toward the bottom. 
[Bright zones represent high average ray density or average squared wave amplitude.] Very impressive and detailed agreement 
is seen comparing the ray density and squared wave amplitude averages. The classical ray density without averaging over initial 
dire ctio n is shown in A in red. Freak events of 4.4cr are recorded in blue in panel D, and 6a events are shown in red. [See 
Fig.~ 



10 for the relevant statistical information. 



Whereas the region shown in Fig. [8] is entirely within the random potential field. Fig. |9] shows wave simulation 
results "before, during, and after" random scattering. Waves are again launched from the top toward the bottom. 
All four panels show the mean squared wave amplitude in greyscale, with panel A showing it unadorned. Panel B 
superimposes the potential field, which is random only in the central region and zero elsewhere. Panels C and D 
superimpose freak event information, with C encoding all 6ct events, and D all 4.4f7 events that occurred during the 
run. Note that very few 4.4cr events, and no 6cr events, occur before the random potential is encountered. Significant 
numbers of these events are found "downstream" of the random potential, but the largest concentration of extreme 
events is located within the refracting region, and especially near the first zone of (smoothed) caustics. Referring back 
to (|9|, we note that the unlucky ship that finds herself in one of the bright lumps in Fig. |9] (where intensity is as high 
as four times the mean) will have approximately 1500 times the probability of encountering a 4.4ct freak wave than 
in a zone of average energy density. A fearsome Scr event is 12,000 times more likely there than in a zone of average 
energy density. 

Figure [To] strongly supports our idea of locally Gaussian statistics averaged over the mean wave density distribution. 
The density for this run (7 = 3.4) is shown in greyscale in the right panel, which also indicates the three regions c, b, 
a over which statistics are collected. Region c is located in the entrance zone before the onset of the random potential; 
region b is located inside the refraction zone and contains the first (smoothed) caustics; and finally region a shows the 
wave behavior after refraction. Naturally, variations are seen from run to run and are largest for the rare events. We 
never obtain sufficient numbers of extreme freak wave (crest height > 6cr) prior to refraction (region c). However we 
see many 6a events during and after refraction. In each case, the dashed line gives the Rayleigh distribution, based 
on the average SWH; the solid red line shows the numerical crest height data, and the dotted line is theory based 
on (10) combined with measurement of the spatial mean density variation I{x,y). [That is, V{I) is determined from 
the mean density data I{x,y), and used in the integral (10 1 to obtain the predicted P{H).] In panel a, the energy 



density distribution V{I) is quite close to Gaussian. A Gaussian distribution of energy density is also a reasonable fit 
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FIG. 9: All the results shown here are obtained from a wave simulation, at 7 ~ 2. 700 randomly chosen plane waves, with 
a range of propagation directions and wavelengths, are superposed and launched from the top toward the bottom. All four 
panels show the mean squared wave amplitude (bright zones indicating high mean squared amplitude), with panel A showing 
this quantity unadorned. Panel B superimposes the potential field, which is random only in the central region. Panels C and 
D superimpose freak event information, with C encoding the 6a events, and D the 4.4(t events that occurred during the run. 



to region b, which is relevant to the derivation of ( 14 1 in the previous Section but not needed here. Note the excellent 



agreement with the Rayleigh distribution in the "undisturbed" region, c. Note too the excellent agreement between 
the data and the predicted curves in all three regions. It is important to note that refraction effects are evident only 
in the tail: the probability distribution of "typical" crest heights {H < 3a) is consistent with Rayleigh in all three 
regions. Consequently, the measured SWH differs in the three regions a, b, and c, by less than 3%. 



Table IV presents some typical data obtained from the linear Schrodinger propagation. For the three typical runs 
shown here, we report AAa, 5a, and 6a and larger events in regions a, b, and c. See Fig.[TO]for an image of the three 
zones. 



V. QUALITATIVE OBSERVATIONS 



We remark briefly on some qualitative and perhaps speculative aspects of the results. Anecdotal reports from 
mariners implicate at least three types of freak waves: 1) The "three sisters", i.e., three large waves in a row heading 
in approximately the mean wave direction, 2) "out of nowhere" waves, which attack suddenly from an unexpected 
quarter, perhaps 45 degrees from the mean wave direction, and 3) the infamous "hole in the sea followed by a wall of 
water". This third kind of wave is reported to be persistent, and perhaps very broad. Certainly, nonlinear effects are 
important in forming and sustaining the third type of wave; it is possible that the linear effects considered here could 
help trigger it. 
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FIG. 10: Log of the crest height probabihty distribution by region (compared with Rayleigh and theory), for 7 = 3.4. The 
dashed line is the Rayleigh distribution, based on the average SWH. The solid line shows numerical data from wave propagation. 
The dotted lines represent the theory based on ( |10| l, together with measurement of the mean density I{x, y). Note the excellent 
agreement with the Rayleigh distribution in the "undisturbed" region, c. 
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TABLE I: For three typical runs shown here, we report events of size at least 4.4(7, 5(7, and 6a in regions a, b, and c. See 
Fig. [To] for an image of the three zones. Runs 1 and 3 have the refracting zone extending throughout regions b and a. Run 
2 has the refracting zone only in region b (as in Fig. [9|. Re gion c is the non-refracting zone preceding first entry into the 
refraction region. Here "ncr pred" is the predicted number of na events (n = 6, 5, 4.4), divided by the expected number based 
on Gaussian statistics over the whole region. Similarly "ncr found" is the theoretically predicted value for this ratio. The 
asterisk indicates that no 6(7 events were seen in the pre-refracting region for the entire run. For a freak index 7 = 1, we see a 
factor of ~ 50 increase in the number of 6a and larger events, a factor of ~ 10 increase in 5a and larger events, and a ~ 4 fold 
increase in 4.4(7 and larger events, compared with the Longuet-Higgins Gaussian seas model in the refracting region b. Run 
1 for 7 — 0.67 {AO ~ 9 degree spread of incident wave directions, and 56^ « 6 degree mean deflection) shows correspondingly 
lower enhancement of freak events over the Gaussian expectations. Run 3, for a freak index of 1.35, shows up to 350 fold 
enhancement of 6a and larger events. 



There is modest evidence for the first two types of wave in our simulations. This is not to say that nonhnearities 
play no role, indeed they must for a complete description. But we find in the simulations that the freak events often 
appear as three sisters, or perhaps five with the middle three being tallest, for the parameters we used. Furthermore, 
we do see rays propagating at high angle to the mean direction; these leave tracks such as the small, fairly bright 
branches seen in Fig. [3] even in the case of a 25 degree spread in the incident waveset. Such streaks are also seen in the 
Schrodinger simulations in Figs. [8] and |9] Waves are seen traveling in the direction of the streaks in the simulations. 
These streaks are typically guided by remnant fold caustics that have survived the averaging and received several 
chance deflections in one direction. They are compact and move at a high angle relative to the mean wave direction. 
We call these streaks and the waves that populate them "runners" . 
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We should also mention another mechanism for lumping, namely the collision or overlapping of two smoothed 
V-shaped fold caustics that form immediately following a cusp. These resemble "rooster tails" seen in the wake of 
power boats, and we call them by that name. Several are seen in Figs. [5] and |9] They cause a sudden increase of local 
average wave action. 

Certainly, the qualitatively reasonable idea that the ray density fluctuations are washed out due to chaotic expo- 
nential instability of the rays is incorrect. What is perhaps even more surprising is that instead of smoothing out, the 
ray density develops finer structure after propagating several times the distance to the first focal regions. 

VI. CONCLUSIONS AND OUTLOOK 

We have seen that refraction of a stochastic Gaussian sea by random eddy currents creates a lumpy spatial energy 
distribution, causing deviations from the initial Rayleigh crest height (or wave height) distribution expected in a 
random wave model. Significant energy lumps survive averaging over wave direction and wavelength, despite the 
chaoticity displayed by individual ray trajectories and the washing out of singularities in the ray dynamics. These 
lumps dramatically increase the probability of freak wave formation, even though parameters associated with low-order 
moments of the wave height distribution, such as the significant wave height or the kurtosis, are almost unchanged. A 
single dimensionless parameter called the freak index 7, defined as the ratio of typical angular deflection in one focal 
distance to the initial angular uncertainty of the incoming waveset, is sufficient to determine the size of the tail in 
the crest height distribution. Although freak wave probability increases with increasing 7, very signiflcant effects are 
obtained already when 7 takes modest values realizable in nature. The increase in extreme freak waves (those of wave 
height greater than three times the significant wave height) is especially spectacular. The number of such waves may 
increase by two or more orders of magnitude when a well-collimated sea (angular spread A6' < 15°) encounters a field 
of strong eddy currents (rms current speed uq = 0.5 m/s). For physically reasonable parameters, good agreement is 
obtained between analytic results based on a small— 7 approximation and numerical simulations. 

In this paper we use statistics averaged over a large area, including all the lumps and streaks, to define the SWH 
and the global statistics. For any reasonable freak index, the low moments of the global distribution, including the 
kurtosis, are scarcely affected, yet as we have shown the far tails of the distribution, i.e., the "freak" events, can 
be greatly enhanced. One could perhaps argue that inside a spatially extended lump one should re-define the SWH 
to adjust to the local energy density and then inside that region the statistics would be the "expected" Gaussian 
with a larger SWH than say some kilometers away. The largest of the lumps are on the order of the correlation 
length of the eddies, which could be thought of as either large or small depending on the wavelength of the seaway. 
There are also much smaller lumps that appear when two streaks cross for example, as well as the fine structure 
mentioned in the previous Section. The basic idea of the non-uniform sampling theory used here is that the seas suffer 
a sudden increase in local wave action when encountering a lump, and they have neither the time nor space to fully 
accommodate this increase through the normal nonlinear evolution. The seas inside a lump have the same wavelength 
but larger amplitude and are therefore steeper, greatly enhancing the probability of a freak wave appearing, just as 
in other scenarios where wave steepness increases. By combining refraction and random wave statistics we retain a 
statistical model for freak wave formation. 

As emphasized in the introduction, the linear model considered here may be regarded a starting point for a more 
sophisticated nonlinear analysis. To leading order in the wave steepness, nonlinear effects may be taken into account 
by replacing the Rayleigh distribution of crest heights in ([S]) with the Tayfun distribution [50]. At this order, the 
only effect of nonlinearity is to create an asymmetry between wave crests and wave troughs, while the distribution of 
wave heights remains unchanged. Numerical evidence suggests that the Tayfun approximation gives accurate results 
for crest height statistics in two dimensions except in the case of a narrow directional spread [21] . 

At higher order, nonlinearity of the wave equation results in exponential instabilities, such as the Benjamin-Feir 
instability |22j for an initial plane wave evolved under the nonlinear Schrodinger equation (NLS). The likely effect 
of such instabilities is an enhancement in the probability of freak wave formation on short to moderate distance 
scales [5j [6] (the Benjamin-Feir time scale is of order {kH)~^ wave periods, where kH is the steepness). On longer 
time scales, nonlinearity should have the opposite effect of reducing wave steepness by transferring energy in the 
hot spots to longer wavelength modes, and resulting in a new equilibrium distribution with a larger significant wave 
height. We do not expect such re-equilibration to be effective when the energy density is varying significantly over 
scales as short as 10 wavelengths. However, further investigation employing a nonlinear model such as a four-wave 
approximation [23 is needed to determine whether refraction and nonlinearity acting together produce more freak 
waves than does either effect separately. Interestingly, recent numerical explorations of nonlinear effects on freak 
wave formation in two dimensions (using nonlinear equations of the Dysthe type) show that these effects are strongly 
dependent on directional spread |211 124] . Nonlinear enhancement of freak wave formation is maximized for long- 
crested waves, which is the same limit in which refraction effects are greatest. Of course, nonlinear effects additionally 



17 



scale with the initial wave steepness (e.g., with the a parameter in the JONSWAP spectrum), while refraction effects 
depend on scattering angle, allowing for nontrivial interplay between the two mechanisms. A number of other issues 
suggest themselves but are not addressed here, such as experimental detection of energy lumpiness, and the rate of 
movement of lumps (due to changing eddy positions and velocities) . 

The essence of this paper is a synthesis of Longuet-Higgins Gaussian seas pj and the refraction model of White 
and Fornberg [9]. From the perspective of wave propagation, it is still a linear theory, and indeed we have checked it 
successfully with a linear Schrodinger propagator. From the point of view of ray dynamics it is decidedly nonlinear 
and even chaotic in the usual sense of exponential sensitivity to small changes in initial data. If the Longuet-Higgins 
Gaussian seas model is "too cold" (too few extreme events are predicted to occur), and the White and Fornberg model 
is "too hot" (freak waves appear periodically at every focal point), then the present combination of Gaussian seas 
and refractive effects may be "just right" . 

We mention that random refraction due to many small eddies is not essential to the main idea presented here, 
which is the notion of random wave statistics in the presence of energy lumps and streaks. Such variations in energy 
or action could have many causes, including a single large eddy or indeed a variable shallow bottom. In fact the latter 
seems the best hope for a real world laboratory to test the ideas in this paper. For a sea incident on a continental shelf 
with bottom depth variations, one might hope to compare the statistics of the incoming wave with the statistics at 
various locations that are predicted to be action maxima or minima. This could hopefully become a semiquantitative 
test of the ideas presented here, i.e., a deviation from Gaussian statistics in the far tails of the distribution. Indeed 
the Canadian waters off Vancouver Island and the Queen Charlotte Islands, an area known for its unruly seas, may 
be one such "laboratory" . 

Finally we express our conviction that nonlinear wave effects are very important, perhaps to every freak wave 
event. However, we believe the sudden buildup of energy or action in the lumps that we have shown exist may be an 
important triggering mechanism for nonlinear evolution. 
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